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Weather  forecasting  in  the  ncsoseale  range  often  remains  problem.it io 
due  to  the  rapid  changes.  With  the  advance  of  high-resolution  imaging 
techniques,  satellite  pictures  taken  in  short  time  intervals  can  he  em¬ 
ployed  to  Improve  interpretation  and  prediction  of  atmospheric  situations 
hv  analyzing  cloud  motion  and  size  parameters.  Knowledge  of  these  para¬ 
meters  allows  estimation  of  wind  fields  as  well  as  the  establishment  ol 
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time  series  in  order  to  model  the  cloud  development  and  extrapolate  from 
past  observations. 

Our  objective  is  to  find  estimators  that  reveal  not  only  information 
on  lateral  (x-y)  displacement,  but  also  on  rotation  of  clouds.  These 
estimators  should  perform  well  for  a  variety  of  meteorological  conditions. 

Currently,  a  variety  of  cloud  displacement  estimators  are  competing 
for  computational  convenience  and  versatility  in  application  [1-63.  Our 
work  concentrates  on  methods  applying  the  Inner  Product  (popularly  called 
cross-correlation)  since  we  hope  to  combine  a  mathematically  elegant 
theory  with  efficient  strategies  for  application.  __ 

V 

Particular  emphasis  will  be  directed  to: 

(1)  Improvement  of  peak  detection  on  the  Inner  Product  surface  by  pre¬ 
filtering  of  the  pictures.  This  method  will  be  called  Accentuated  Corre¬ 
lation  (AC).  An  analysis  of  the  given  signals  and  their  properties  in  the 
spatial  and  frequency  domain  leads  to  the  theoretical  background  of  AC. 

AC  will  adopt  ideas  developed  for  the  one-dimensional  case  [7]  and  incor¬ 
porate  a  recently  suggested  exponential  type  of  pre-filtering  [8].  It  can 
be  shown  that  pre-filtering  with  a  truncated  inverse  filter  leads  to 
better  results  than  the  exponential  type. 

(?)  Employing  Inner  Product  methods  on  the  polar-coordinate  contour  func¬ 
tions  of  two  clouds  to  determine  their  relative  rotation. 

At  the  time  of  the  workshop,  results  of  the  suggested  methodr  will 
be  presented  for  the  case  of  simple,  single  clouds.  An  extension  for 
more  complex  cloud  situations  is  currently  under  work. 


Accession  Far 

wris  «A*i 

DDC  TAB 

Unannounced 

Justification 


:  By _ _ 

Distribution/ 

!  Availability  Codes 
Avail  and/or 
iDlst  special 


Unclassified 


MCuetTv  cl  ami  pic  axiom  or  >»"  rAoenrhm  om  Ent»r*o 


1)  Introduction 

Weather  forecasting  in  the  mesoscale  range  often  remains  problematic 
due  to  the  rapid  development  of  disturbances.  With  the  advance  of  high- 
resolution  imaging  techniques,  satellite  pictures  taken  in  short  time 
intervals  can  be  employed  to  analyze  motion  and  development  of  clouds. 
Knowledge  of  these  parameters  allows  for  the  estimation  of  windfields  as 
well  as  for  the  prediction  of  certain  types  of  severe  weather  situations. 

Our  objective  is  to  find  estimations  that  reveal  information  not 
only  of  lateral  (x-y)  displacement,  but  also  of  cloud  rotation  and  size 
changes.  The  scope  of  application  might  well  be  extended  to  images  of 
moving  objects  in  general. 

Currently  a  variety  of  lateral  cloud  displacement  estimators  are 
competing  for  computational  convenience  and  versatility  in  application: 
the  film  loop  technique  [1],  the  "correlation"  method  [2,3],  binary 
matching  [2],  clustering  analysis  [4,5]  and  relaxation  techniques  [6]. 

Our  work  concentrates  on  methods  applying  the  inner  product  (often 
misinterpreted  as  cross-correlation) ,  since  we  hope  to  combine  a  mathe¬ 
matically  elegant  theory  with  efficient  strategies  for  application  in 
not-too-complex  cloud  situations. 

In  this  paper,  we  shall  introduce  a  new  pre-filtering  technique 
and  suggest  a  method  to  identify  the  rotation  of  the  cloud.  A  rotating 
cloud  very  often  indicates  a  severe  weather  situation.  Conventional 
inner  product  methods  suffer  from  poor  maximum  peak  resolution  due  to 
cloud  shape  changes,  rotation  and  other  disturbances  within  the  frame 
of  interest.  To  improve  maximum  peak  detection,  pre-filtering  of  the 
input  images  can  be  performed  prior  to  the  inner  product  operation. 


For  the  case  of  one-dimensional  signals,  it  has  been  shown  [7]  that  inverse 
filtering  based  on  the  knowledge  of  the  signal  spectrum  leads  to  improve¬ 
ment  in  certain  instances.  Recently  a  complex  exponentiated  type  of  pre¬ 
filtering  was  suggested  [8],  The  method  we  introduce  will  be  called 
accentuated  inner  product  (AIP) .  The  pre-filter  is  a  two-dimensional 
quasi-difference  filter  with  its  coefficients  determined  from  the  ’auto- 
inner-  product'  matrix  by  solving  the  Yule-Walker-equations .  In  comparison 
with  the  unfiltered  and  the  complex  exponentiated  inner  product,  our  em¬ 
pirical  studies  show  the  superiority  of  the  AIP. 

In  estimating  the  rotation  of  a  cloud,  we  introduce  a  method  based 
on  the  cross-correlation  estimate  of  the  polar-coordinate  envelope 
functions. 

Finally  we  suggest  the  Incorporation  of  simple  motion  models  to 
derive  initial  guesses  prior  to  the  analysis  of  the  currently  received 
satellite  picture  in  order  to  cut  computation  time  and  allow  for  fast 
identification  of  severe  storm  situations. 

2)  Signal  Analysis  and  the  Inner  Product 

As  conventional  prerequisites,  the  given  images  are  assumed  to  be 
discrete  in  a  regular  cartesian  grid,  and  the  pixels  containing  a  dig¬ 
ital  brightness  value,  for  example  any  integer  number  between  0  (black) 
and  255  (white).  For  displacement  analysis,  all  image  frames  have  to 
be  navigated  to  some  reference  mark,  and,  if  possible,  the  brightness 
of  the  clouds  adjusted  for  sun  angle  changes.  With  the  exception  of 
snow-covered  terrain,  the  background  usually  is  darker  than  the  clouds. 

It  is  removed  by  subtracting  the  threshold  value  from  the  picture  values. 
Furthermore,  we  assume  the  cloud  of  interest  is  completely  defined  within 


the  image  frame. 
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Since  the  cloud  image  sequence  {x(n,m)>  then  is  zero  outside  some 
interval  0  <  n  <  N-l  and  0  <  m  <  M-l,  it  belongs  to  the  class  of  "finite 
area  sequences".  Finite  area  sequences  are  energy  signals,  i.e. 

0  <  | |x(n,m) | | 2  <  •  (1) 

with  the  energy  of  (x(n,m)}  defined  as 

2 

|  |x(n,m)  |  |  -<x(n,m),  x(n,m)> 

00  00 

-  Z  Z  x2(n,m) .  (2) 

H«-CD  IQS— 00 

These  signals  describe  the  space  with  the  inner  product  of  two 
sequences  {x(n,m)}  and  {y(n,m)},  defined  as 

K(0,0)  -  <x(n,m) ,  y(n,m)> 

M-l  N-l 

*  Z  Z  x(n,m)y (n,m) .  (3) 

m»0  n-0 

If  the  sequence  {y(n,m)}  is  displaced  by  v  and  p,  we  obtain 

M-l  N-l 

K(v,p)  -  Z  Z  x(n,m)y(n-v,m-p) .  (4) 

m*0  n*0 

Note  that  the  cross-correlation  is  defined  as  a  statistical  expected 
value  type  of  operation  involving  the  Joint  probability  density 
function.  Energy  signals  are  neither  stationary  nor  ergodic.  Thus, 
we  shall  refer  to  displacement  estimators  involving  functions  like 
eq.  (4)  as  the  inner  product  method.  If  v,p  in  eq.  (4)  are  varied 
from  0  ...  N-l,  0  ...  M-l,  there  is  one  set  of  indices  v,u  such  that 


K(*>  is  maximum. 
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For  the  case  that  {y(n,m)>  is  a  shifted  version  of  {x(n,m)},  i.e. 

{y(n,m) }  -  {x(n-n0,m-ra0) } ,  (5) 

K(*)  is  maximum  for  v  =  n^,  p  *  m^.  For  finite  area  sequences,  the 
Discrete  Fourier  Transform  can  be  defined.  With  (X(n,m)}  and  {Y(n,m)} 
the  DFT  sequences  of  {x(n,m)}  and  {y(n,m)},  Parseval’s  theorem 


.  1  N~1  MI1  v/  .  i2ir(um+vn)/MN 

K(v,p)  *  -ng  E  E  X(n,m)Y(n,m)e 

n=0  m*0 


(6) 


shows  that  displacement  estimation  can  also  be  performed  by  estimating 
the  phase  of  the  cross  spectrum.  However,  due  to  the  current  lack  of 
suitable  unwarping  algorithms,  the  evaluation  of  (4)  is  preferred  to 
the  spectral  analysis. 

In  Fig.  1,  the  inner  product  of  a  one-dimensional  square  wave 
signal  (x(k) }  and  its  shifted  version  (x(k-nQ)}  is  shown,  both  some¬ 
what  resembling  a  cloud  profile.  Note  that  the  resulting  triangular 
function  has  its  maximum  at  t  •  n^,  but  is  twice  as  wide  as  the  orig¬ 
inal  function.  This  effect  is  well  known  from  the  convolution  pro¬ 
perties.  The  inner  product  is  at  least  as  wide  as  the  widest  of  both 
the  input  functions. 

Since  a  broad  cross  spectrum  is  equivalent  to  a  narrow  inner  product 
peak,  whitening  filters  of  various  kinds  were  suggested  [7]  based  on  the 
knowledge  of  the  signal  spectrum.  If  that  spectrum  is  invertible,  a 
filter  designed  to  approximate  the  inverse  spectrum  would  "whiten"  the 
input  signals  prior  to  the  inner  product  operation.  Recently  a  ronlinear 
complex  exponentiation  pre-filter  was  suggested  [8],  of  the  form 


x(n,m)  exp(-ipx(n,m)) 


(7) 
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with  p  some  constant.  This  pre-filter  results  in  an  absolute  value  inner 
product 
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with 


w  =  p(x(n,m)-y(n+v,m+|j)) . 


The  exponential  pre-filter  was  introduced  by  studying  optical  trans¬ 
parency  properties. 

With  the  series  representation 


e 


-ix 


2  3 

ix  x  lx 

11  ~  21  3!  ***  ’ 


(9) 


the  exponential  pre-filter  creates  a  series  of  harmonics  and  thus  broadens 
the  spectrum.  However,  as  our  results  show,  it  is  rather  sensitive  to 
object  changes. 


3)  Accentuated  Inner  Product 

In  Fig.  2  it  is  shown  that  a  differentiator-type  pre-filter  for  the 
example  of  Fig.  1  leads  to  a  delta-function  at  t  =  ng-  Note  that  the 
frequency  response  of  the  differentiator  is  a  high-pass,  thus  this  filter 
removes  the  dominating  low  frequencies  of  energy  signals  and  serves  as 
a  whitening  filter  to  some  extent. 

We  Introduce  the  two-dimensional  quasi-difference  filter 


1  1 

y(n,m)  -  l  Z  a(i, j)x(n-i,m-j)  (10) 

i-0  j*0 

which,  in  effect,  is  a  two-dimensional  moving  average  filter  of  order  two. 
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Its  coefficients  are  chosen  by  the  predictor  filter  design  technique. 
Thus,  the  two-dimensional  Yule-Walker  equations 


”r(0,0)  r (0,1)  r(l,0)  r(l,l) 

Y(0,0) 

1 

r (0,-1)  r(0,0)  r (1,-1)  r(l,0) 

Y (0,1) 

0 

r (-1,0)  r(-l,l)  r (0,0)  r(0,l) 

Y(1.0) 

0 

r(-l,-l)r(-l,0)  r(0,-l)  r(0,0)_ 

Y(l,l)_ 

_0_ 

with  r(j,k)  =  <x(n,m),  x(n-j,m-k)>  ,  (11) 

are  solved  for  the  vector  y.  Then,  a(j,k)  *  -y(j ,k)/y(0,0)  for  j ,k  +  0, 
and  a(0,0)  =  1.  Note  that  is  the  first  column  of  the  inverted  matrix 
R.  The  design  of  a  prediction  filter  of  higher  order  would  be  detri¬ 
mental  to  the  result,  since  our  given  signals  are  in  general  not  of  the 
auto-regressive  type,  i.e.  have  no  poles.  The  zero  of  the  quasi-dif¬ 
ferentiator  just  cancels  the  high  amount  of  low-frequencies  character¬ 
istic  to  finite  area  sequences. 

Results 

In  Fig.  3  the  inner  product  surfaces  are  shown  for  no  pre-filtering 
(top),  exponentiated  pre-filtering  (middle),  and  the  AIP.  The  plots  in 
the  left  column  were  obtained  by  performing  the  inner  product  of  a  cloud 
and  its  shifted  replica.  The  plots  of  the  right  column  correspond  to 
cloud  pairs  with  actual  changes  as  observed  in  a  3-minute  interval. 

For  the  exponentiated  filter,  the  factor  p  ■  5  seemed  to  give  best  re¬ 
sults.  Although  the  delta  peak  of  the  exponentiated  inner  product  for 
the  unchanged  cloud  is  quite  impressive,  the  AIP  method  seemed  to  be 
superior  for  the  case  of  changing  clouds.  Since  the  output  is  deter¬ 
mined  for  a  regular  grid,  further  improvement  for  estimating  the  dis¬ 
placement  could  be  achieved  by  interpolating  over  the  values  in  the 
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vicinity  of  the  maximum.  A  bicubic  spline  seems  to  be  feasible  to  obtain 
the  maximum  between  grid  points  by  solving  its  derivative  for  zero. 

4)  Rotational  Displacement  Estimation 

As  indicated  above,  often  it  is  desired  to  obtain  an  estimate  of  the 
rotation  of  the  cloud,  be  it  to  improve  the  lateral  displacement  estimate 
by  turning  the  image  to  a  best  match,  or  be  it  for  gaining  information 
from  the  amount  of  rotation  for  weather  analysis.  In  the  following  we 
assume  rotation  around  a  point  inside  the  object,  possibly  the  center  of 
brightness  gravity.  The  detection  of  the  actual  axis  of  rotation,  how¬ 
ever,  still  remains  a  problem,  although  in  many  cases  the  center  of 
gravity  or  some  other  reference  point  commonly  recognized  in  all  images 
of  the  cloud  sequence  is  sufficient.  Fig.  4  shows  a  contour  plot  of  a 
cloud  for  three  different  gray  levels  (L=3,  12  and  35,  after  subtraction 
of  the  background  level) .  From  these  contours  and  the  coordinates  of 
a  given  point  inside,  a  polar-coordinate  envelope  is  determined  by 
obtaining  the  distance  (radius)  to  the  contour  for  an  equidistant  angle 
domain.  Fig.  5  shows  the  three  envelope  functions  for  the  cloud  in 
Fig.  4  over  the  domain  0°-360°,  0°  being  parallel  to  the  x-axis. 

The  envelope  curves  are  periodic  with  periodicity  2 it.  If  a  cloud 
rotates,  the  envelope  functions  shift.  Hence  the  amount  of  rotation  can 
be  estimated  by  estimating  the  cross-correlation  of  the  envelope  func¬ 
tions.  Fig.  6  shows  the  result  of  a  Blackman-Tukey-type  cross-correla¬ 
tion  of  an  L  =  35  envelope  with  its  rotated  successor  indicating  a 
maximum  at  lag  10,  i.e.  20°  rotation  for  each  lag  representing  2°.  Note 
that  the  peak  is  more  discernible  for  the  L  =  12  envelopes.  The  pro¬ 
posed  method  is  currently  tested  for  its  practicability  in  estimating 


rotational  displacement.  However,  the  Fourier  transform  of  the  envelope 
function  might  be  of  great  value  for  identifying  the  same  cloud  in  con¬ 
secutive  images. 


5)  Cloud  Motion  Modeling 

A  problem  common  to  all  inner  product  methods  is  the  large  amount 
of  computing  time  necessary  to  perform  the  required  operations.  Com¬ 
puting  time  can  drastically  be  reduced  if  an  estimate  of  the  cloud 
motion  exists  beforehand.  Then,  the  output  window  can  be  set  as  small 
as  possible  around  the  expected  surface  maximum.  On  our  PDP  11/60,  the 
computing  time  was  reduced  to  less  than  a  minute  for  obtaining  a  6  by  6 
output  surface  from  128  by  64  pictures,  pre-filtering  included. 

Feasible  predictors  are  currently  under  research.  Our  goal  is  to 
implement  some  scalar  cloud  parameters  such  as  area  size,  brightness 
concentration,  etc.  These  are  easily  obtainable  and  should  be  combined 
with  the  displacement  vectors  to  low-order  models  that  allow  for  quasi- 
real-time  cloud  state  identification  and  short-time  prediction. 


t 
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Fig,  3.  Inner  Product  surfaces  for  unfil 
(middle)  and  accentuated  (bottom)  inputs 
for  changed  cloud. 


Fig.  6.  Cross-correlation  estimates  of  envelope  functions. 


